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ABSTRACT 


In  this  thesis,  the  performance  and  accuracy  of  explicit,  semi-implicit,  and 
Hybrid  Eulerian-Lagrangian  Semi-Implicit  (HELSI)  time-integration  methods  for  use 
in  atmospheric  modeling  are  examined.  Four  test  cases  are  analyzed:  A  density 
current,  an  inertial  gravity  wave,  a  rising  thermal  bubble,  and  a  hydrostatic  mountain 
wave.  Strict  attention  is  paid  to  computational  time,  stability  criteria,  and  accuracy. 
The  project  aims  to  show  increased  efficiency  using  the  HELSI  method  over  semi- 
implicit  methods,  which,  in  turn,  should  be  better  than  the  split-explicit  methods 
currently  used  in  mesoscale  models  such  as  WRF,  COAMPS,  and  the  German  LM 
model.  This  increase  in  efficiency  allows  for  valuable  computational  resources  to 
be  used  for  other  purposes,  such  as  improved  data  assimilation,  increased  spatial 
resolution,  or  more  detailed  physics. 
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I.  INTRODUCTION 

A  numerical  weather  prediction  (NWP)  modeling  system  is  a  method  of  pre¬ 
dicting  the  future  state  of  the  atmosphere  using  data  collected  about  the  current  state 
of  the  atmosphere  and  known  information  about  atmospheric  behavior.  These  behav¬ 
ioral  properties  are  approximated  using  a  mathematical  model,  called  the  governing 
equations,  and  by  physical  parametrizations. 

An  operational  NWP  modeling  system  is  a  multi-component  system.  Processes 
for  accurate  data  assimilation,  forecast  integration,  and  physical  parametrizations  are 
vital  to  produce  an  accurate  forecast.  Improving  any  one  of  the  components  should 
increase  the  overall  accuracy  of  the  NWP  system,  received  as  “NWP  guidance”  by 
forecasters,  and  should  thus  help  produce  a  better  operational  forecast. 

Improvements  to  any  component,  however,  generally  require  increased  compu¬ 
tational  expense.  When  improvements  to  NWP  systems  are  made,  timeliness  must 
also  be  maintained.  Increasing  computational  requirements  must  be  matched  by  in¬ 
creased  computational  capacity,  or  increased  computational  efficiency,  or  both.  Lim¬ 
ited  budgetary  resources  demand  that  significant  attention  ought  to  be  placed  on  the 
computational  efficiency  approach.  One  area  where  a  significant  improvement  in  effi¬ 
ciency  may  be  achieved  is  the  numerical  time-integration  methods  used  to  march  the 
governing  equations  forward  in  time,  and  it  is  the  purpose  of  this  thesis  to  advance 
the  prospect  of  implementing  more  efficient  methods  in  Department  of  Defense  NWP 
systems. 

Many  current  non-hydrostatic  mesoscale  NWP  models  such  as  the  Weather 
Research  and  Forecasting  (WRF)  model  and  The  LI.S.  Navy’s  Coupled  Ocean  At¬ 
mosphere  Mesoscale  Prediction  System  (COAMPS)  use  split-explicit  time  integra¬ 
tion  methods.  Semi-implicit  methods  suggested  by  Giraldo  [1]  offer  a  significant 
improvement  in  efficiency.  Hybrid  Eulerian-Lagrangian  Semi-Implicit  (HELSI)  time- 
integrators  may  offer  even  further  improvements. 


1 


Here,  the  accuracy  and  efficiency  of  explicit,  semi-implicit,  and  HELSI  meth¬ 
ods  are  explored.  Four  test  cases  are  examined:  a  rising  thermal  bubble,  a  linear 
hydrostatic  mountain  wave,  a  density  current,  and  an  inertia-gravity  wave.  The 
various  time  integrators  are  used  to  calculate  the  perturbations  in  Exner  pressure, 
velocity,  and  potential  temperature  for  each  of  these  test  cases.  Careful  attention  is 
paid  to  wallclock  computational  time,  stability  criteria,  and  accuracy. 
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II. 


BACKGROUND 


A.  GOVERNING  EQUATIONS 

For  this  study,  equation  set  1  of  Giraldo  and  Restelli  [2]  is  used.  This  equation 
set  is  the  non-conservative  form  of  the  Euler  equations.  It  is  the  two-dimensional 
version  of  the  equations  currently  used  in  operational  NWP  models  such  as  COAMPS 
and  the  German  LM  model. 

1.  Equation  Set 


dn  R 

— — b  u  ■  V7T  -| - 7 rv  •  u 

at  cv 

dii 

— - b  U  ■  V«  +  C„dV 7T 

at 

39  ™ 

at 


0 

— gk  +  //.V2w 


(2.1) 


Here,  7r  =  ^  v  is  the  Exner  pressure,  u  =  ( u ,  iv)T  is  the  2-dimensional  velocity 

held,  and  6  =  ^  is  the  potential  temperature.  The  solution  vector  is  (7 


B.  SPATIAL  DISCRETIZATION 

A  finite  element  spatial  discretization  on  quadrilateral  elements,  as  defined  by 
Giraldo  and  Restelli  [2],  is  presented.  This  decomposition  of  the  global  domain  into 
a  multitude  of  smaller  domains  is  designed  to  fully  exploit  the  multiple  processor 
architecture  currently  used  in  high  performance  computing  [1],  In  order  to  solve  the 
equations  numerically,  they  must  first  be  recast  in  the  form 

I  =  s '<«>  <2-2> 

where  S(q )  contains  the  source  terms  of  the  governing  equations.  Then,  the  global 
domain,  fi,  is  broken  up  into  Ne  non-overlapping  quadrilateral  elements  such  that 

Ne 

Q  =  U  Qe. 

e=l 
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In  order  to  perform  calculus  operations,  a  non-singular  mapping  x  =  'I'  (£) 
from  the  physical  Cartesian  coordinate  system  x  =  (x,  z )  G  f le  to  the  reference 
coordinate  system  £  =  (£,  rj)  which  is  defined  in  each  element  such  that  (£,  rj)  G 
[— 1,+1]2  in  each  element.  Also  used  is  the  transformation  Jacobian,  J  =  which 
is  associated  with  the  local  mapping,  This  mapping  is  then  used  to  define  the 
local  representation  of  q.  the  solution  vector,  and  the  approximation  of  the  calculus 
operations. 

The  structure  of  the  reference  element,  I,  which  is  spanned  by  £  G  [ — 1,  l]2, 
makes  a  natural  representation  of  the  local  element-wise  solution  q  by  an  TVth  order 
polynomial  as  follows: 

K 

qn(x)  =  5 ZMx)qN(xk ) 

k=  1 

where  Xk  represents  K  —  (N  +  l)2  grid  points  and  ipk{x)  reflects  the  associated 
multivariate  Lagrange  polynomials;  these  polynomials  are  the  basis  functions  used  in 
standard  finite  elements.  The  square  structure  of  the  reference  element  allows  the 
representation  of  the  Lagrange  polynomial  by  a  tensor-product 


^ k(x )  =  hi{£{x))hj{r]{x)),  (2.3) 

where  i .  j  =  0, ...,  N. 

For  the  grid  points  (£i,rjj)  the  Legendre-Gauss-Lobatto  (LGL)  points  are  cho¬ 
sen,  given  as  the  tensor-product  of  the  roots  of 

(i  -  emo  =  o, 


where  Pn(0  is  the  Nth  order  Legendre  polynomial.  This  will  simplify  the  algorithm, 
as  the  LGL  points  will  be  used  elsewhere. 

The  one-dimensional  Lagrange  polynomials,  /q(£)  are 


hi(0 


i  (i-aa) 

N(N  +  1)  (£ -6)^(6)’ 
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and  likewise  for  hj{rj)  [3]. 

Having  chosen  the  LGL  points  makes  the  approximation  of  local  element  in¬ 
tegrals  straightforward 


where  J  represents  the  local  Jacobian  of  the  transformation  between  fle  and  I,  and 

c j(£i)  and  uj(rjj)  are  the  Gaussian  quadrature  weights, 

2  /  1 
=  N(N+  1) 

associated  with  the  one-dimensional  LGL  quadrature. 

In  the  spectral  element  (SE)  method,  which  is  a  high-order  finite  element 
method,  a  polynomial  expansion 


K 

QN{x,t)  =  Y,^k{x)qk(t)  (2.4) 

k= 1 

is  used  to  approximate  q.  The  variational  statement  of  equation  2.2  is:  find  q  G 


Hl{yL)  V  i/j  G  H1  such  that 


(2.5) 


where  H 1  (fl)  is  dehned  as  the  space  of  all  C°  continuous  functions  with  functions  and 
first  derivatives  belonging  to  L2(fl),  the  square  integrable  functions,  i.e. 


J  \q\2dQ  <  oo,\/q  G  fh 

Simplifying  Eq.  (2.5)  by  virtue  of  the  SE  discretization  yields  the  global  matrix 
problem 


Ottt  t  ,  ,_i  ^  R  , , 

—  +  ujMj  1Dijttj  +  —njMj  1DjJuj 

du T  t  i  i 

—  +  XljMj  1DIjUj  +  CpdIMI  1DIj1Tj 

d^  +  ujMr1DIJeJ 


0 


-gk  -  \iLuUj 


— h LijOj 


(2.6) 
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where 


Ne 

M  —  f\  Me 

e=l 

is  the  global  mass  matrix 


Ne 

D  =  A  De 


e=l 


is  the  global  differentiation  matrix  and 


Ne  Ne 

L=  f\Le=  f\  /  V^iV^-dfie 

e=l  e=l 


is  the  global  Laplacian  matrix.  In  the  SE  method,  the  local  element-wise  matrices 
are  built  Erst, 

M'  =  (2.7) 

Dtj  =  <|  (2.8) 

£'■  =  E  <  I  I  V*(*i)  •  V^(*(),  (2.9) 

i=l 

where  zce  are  the  quadrature  weights,  Je  is  the  Jacobian,  i  =  1, ...  ,K  are  the  local 
element  grid  points,  e  —  1, Ne  are  the  elements  covering  the  global  domain,  and 
I,  J  =  1, . . . ,  iVp  are  the  global  grid  points.  Finally,  the  direct  stiffness  summation 
(DSS)  operator 

Ne 

A 

e=l 

is  used  to  map  the  local  element  grid  points  i  =  1, . . . ,  K  and  elements  e  =  1 , ,Ne 
to  the  corresponding  global  grid  points  via  the  mapping  (i,  e)  — >  (/)  where  /  = 
1, ...  ,Np  are  the  global  grid  points. 


C.  TIME-INTEGRATION 

A  numerical  time-integrator  is  a  numerical  method  used  to  approximate  the 
time-derivative  of  the  equation 

=  S(q).  (2.10) 
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Here,  q  =  (n,  ur ,  9)T  is  the  solution  vector  of  Equation  2.1,  and  S(q)  is  the  spatial 
discretization  of  the  right  hand  side  of  Equation  2.1. 

There  are  two  basic  types  of  time- integrators:  Eulerian  and  Lagrangian.  [4] 

1.  Eulerian  Methods 

Eulerian  methods  use  a  fixed  frame  of  reference.  As  an  example,  equation  2.10 
is  written  in  Eulerian  form.  As  another  example,  the  advection  equation  in  Eulerian 
form  is 


-^  =  -n-Vq  (2.11) 

where  u  is  the  velocity  vector 

2.  Lagrangian  Methods 

In  contrast,  Lagrangian  methods  use  a  moving  frame  of  reference.  Therefore, 
equation  2.10  written  in  Lagrangian  form  is 

I  =  s ‘<«>  <2-12> 

where 


and 


d 

dt 


d 

dt 


+  u  ■  V. 


u  = 


dx 

dt 


The  advection  equation  written  in  Lagrangian  form  can  be  simplified  to 


^=0. 

dt 


(2.13) 


(2.14) 


(2.15) 


3.  Courant-Friedrichs-Levy  Condition 

A  vital  component  to  a  discussion  on  time-integrators  is  the  Courant  number, 
which  measures  the  amount  of  information  traversing  a  grid  cell  (Ax)  in  a  given 
timestep  (At). 
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For  example,  the  one-dimensional  Courant  number  is  defined  to  be 


C  = 


uAt 

Ax  ’ 


(2.16) 


The  Courant-Friedrichs-Levy  (CFL)  condition  is  the  largest  Courant  number 
that  a  time-integrator  can  use.  In  general,  increasing  the  Courant  number  is  a  central 
goal  to  increasing  the  overall  efficiency  of  the  time-integration  method. 


III. 


TIME-INTEGRATION  METHODS 


A.  EXPLICIT  TIME-INTEGRATORS 

Explicit  Time-Integrators  are,  by  far,  the  simplest  of  the  methods  examined. 
They  have  been  proven  to  be  very  accurate.  They  are,  however,  quite  inefficient.  One 
particularly  nice  feature  is  that  only  knowledge  of  the  solution  at  previous  time  steps 
is  needed.  In  general,  an  explicit  time-integrator  is  of  the  form 

Me  Me 

<f+1  =  +  7A  (EMr).  (3.1) 

2=0  2=0 

1.  Backwards  Difference  Formula 

The  second  order  backwards  difference  formula  (BDF2)  is 

9“+1  =  \qn  +  Iq"-1  +  \u(2S(qn)  +  S(qn~1)).  (3.2) 

2.  Leapfrog 

The  second  order  leapfrog  formula  (LF2)  is 

qn+i  =  qn- 1  +  2At(S{qn)).  (3.3) 

3.  Leapfrog  Stability  Analysis 

Due  to  the  existence  of  Leapfrog’s  computational  mode,  a  Robert- Asselin  time- 
filter  is  necessary  to  preserve  stability.  At  each  timestep,  the  following  filter  is  applied 


qn  =  qn  +  e(qn+1  -  2 qn  +  qn~l)  (3.4) 

where  q  is  the  time- filtered  variable  and  e  is  the  time-filter  weight.  A  brief  linear 
stability  analysis  of  a  simple  equation  demonstrates  this  necessity. 


a.  No  Time-filter 

Consider  the  equation 


dq 

~dt 


=  ikq 
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(3.5) 


whose  solution  is  known  to  be 


q  =  Aeikt  (3.6) 

where  A  is  any  constant.  If  the  equations  are  discretized,  then,  at  timestep  n,  t  =  nAt, 

qn  =  AeiknAt  =  A(eikAt)n  (3.7) 

and 

S(q)n  =  ikqn  =  ikA(elkAt)n.  (3.8) 

Consider  the  case  where  there  is  no  time-filter,  ie  e  =  0,  then 

A(eikAt)n+1  =  A(eikAt)n~1  +  2A  tikA(eikM)n  (3.9) 

or 

An+1  =  A”-1  +  2ikAt\n  (3.10) 

where  A  =  elkAt  is  the  amplification  factor.  If  |  A  |  >  1,  the  solution  grows  exponentially 
fast,  and,  thus,  is  unstable.  Simplifying  yields 

A2  -  2ikAt\  -1  =  0  (3.11) 

which  has  the  solution 

A  =  ikAt  ±  yJ-(kAt)2  +  1  (3.12) 

Note  that  the  amplification  factor  can  take  on  two  values:  A  =  ikAt  + 
^J-(kAt)2  +  1,  which  is  the  physical  solution,  while  A  =  ikAt  —  ^J-(kAt)2  +  1  is 
called  the  computational  mode.  In  order  to  examine  |A|,  two  cases  are  considered: 

|  A; At  |  >  1  and  |A;At|  <  1 
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If  |  A; At |  >  1  the  amplification  factor  is  purely  imaginary.  Thus,  the 
magnitude  of  the  two  modes  are  different.  The  magnitude  of  the  amplification  factor 
of  the  physical  solution  is  ^2( [kAt +  2kAt^(kAt)2  —  1  —  1,  which  is  always  greater 
than  one.  Therefore,  the  method  is  always  unstable  here.  For  completeness,  the 
magnitude  of  the  amplitude  factor  of  the  computational  mode  is 
\j2(kAt)2  —  2kAt\J  (kAt)2  —  1  —  1. 

When  \kAt\  <  1,  the  term  inside  the  square  root  is  real  while  the  other 
term  is  imaginary,  and  the  magnitude  of  the  amplification  factor  for  both  the  physical 
solution  and  computational  mode  is 

(kAt)2  -  (kAt)2  +  1  =  1  (3.13) 

which  is  neutrally  stable. 

Since  the  computational  mode  is  as  large  as  the  physical  solution,  non¬ 
linear  interactions  can  affect  the  physical  solution  and  cause  instabilities  [4],  This 
problem  is  corrected  with  time-filtering.  It  is  necessary  to  note,  however,  that  this 
will  reduce  the  accuracy  of  leapfrog  to  first  order. 

b.  Time- filtered 

Now  consider  the  case  where  the  time-filter  is  applied.  Starting  with 
equation  3.3  and  3.4,  the  definition  of  the  time-filtered  leapfrog  method  becomes 

qn+ 1  =  ~qn~l  +  2A  t(S(qn)).  (3.14) 

Using  the  same  equation 

=  ikq  (3.15) 

gives 

q  =  Aeikt  (3.16) 
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for  the  non-filtered  solution  and 


q  =  Aeikt  (3.17) 

for  the  time-filtered  solution.  Substituting  A  =  elki .  and  dividing  by  An_1,  the  time- 
hltered  leapfrog  method  becomes 

AX2  =  A  +  2AikAt\  (3.18) 

where 

A  =  A  +  e(AX  -  2A  +  AA"1)  (3.19) 

which  follows  directly  from  the  definition  of  the  time-filter.  Simplifying  and  solving 
using  the  quadratic  equation  yields 


A  =  e  +  A tik  ±  \J (e  +  A tik)2  +  1  —  2e  —  2A tike.  (3.20) 

For  a  more  rigorous  treatment,  see  Appendix  B.  The  magnitude  of  both 
the  phsical  solution  and  computational  mode  are  plotted  using  a  value  of  e  =  .05. 
Clearly,  the  computational  mode  is  damped,  and,  thus,  the  problems  of  nonlinear 
interactions  have  been  addressed. 
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a)  b) 


Figure  1.  The  stability  of  the  explicit  leapfrog  time-integrator.  Figure  a)  has  no 
time-filter,  while  figure  b)  has  a  time-filter  weight  of  e=.05.  The  solid  lines  represent 
the  physical  solutions  while  the  dashed  lines  represent  the  computational  modes. 


4.  Runge-Kutta  Methods 

The  Runge-Kutta  (R-K)  methods  examined  are  part  of  a  class  of  methods 
known  as  strongly  stability  preserving  (SSP)  methods  proposed  by  Cockburn  and 
Shu  [5]  and  Ruuth  and  Spiteri  [6].  They  are  third  order,  multi-stage  methods.  These 
methods  offer  larger  stability  regions  than  the  previous  explicit  methods.  The  trade 
off,  however,  is  that  more  computations  are  required  per  timestep.  The  first  of  these 
methods  is  the  three-stage  third  order  R-K  method  (RK3) 


qW  =  qn  +  AtS(qn )  (3.21) 

qm  =  +  1«(1)  +  jAtS(«<») 

«(n+1)  =  \qn  +  ^l2)  + 


13 


The  second  is  the  four-stage  third  order  R-K  method  (RK34) 


qa)  =  <j"+AtS(q")  (3.22) 

«(2>  = 

9,3)  =  +  K(2)  +  AiS(,(21) 

3  3  o 

q(n+i)  =  q(3)  +  ^Af5(q(3)). 

The  final  is  the  five-stage  third  order  R-K  method  (RK35) 

qlx)  =  qn  +  aiAtS(qn)  (3.23) 

qM  =  qW  +  a2AtS(qW) 

q(3)  =  a3qn  +  a4q(2)  +  a5AtS(q(2) ) 

g(3)  =  a6q”  +  a7q(3)  +  a8Af5(q(3)) 

g(n+i)  _  0gg(4)  _|_  a10q(2)  +  anAf5(q(4)) 
whose  coefficients  a*  are  listed  in  Appendix  A. 

B.  SEMI-IMPLICIT  TIME-INTEGRATORS 

Implicit  time-integrators  are  more  complicated  than  explicit  integrators  be¬ 
cause  they  require  knowledge  of  the  solution  at  the  current  timestep  as  well  as  the 
previous  ones.  The  benefit  is  that  they  allow  for  a  higher  Courant  number,  espe¬ 
cially  important  in  the  calculation  of  fast  propagating  waves.  In  general,  an  implicit 
time- integrator  is  of  the  form 

M/  Mi 

qn+1  =  E  +  7  At  E  PiS(qn~l)  (3-24) 

i= 0  i=— 1 

The  semi-implicit  (SI)  method  treats  some  terms  implicity  and  others  explic¬ 
itly.  The  SI  method  of  Giraldo  is  presented  [1]. 
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A  SI  time-integrator  splits  the  source  terms  into  linear  and  nonlinear  terms 


=  { S(q )  -  SL(q)}  +  [SL(q)]  (3.25) 

where  the  terms  in  braces  are  time-integrated  explicitly,  while  the  terms  in  the  square 
brackets  are  time-integrated  implicitly.  L  represents  the  linearization  of  S.  Therefore, 
nonlinear  terms  are  integrated  explicitly  and  the  linear  terms  implicitly.  Finally,  5  =  0 
for  a  fully  explicit  method  and  6  =  1  for  a  semi-implicit  method.  For  this  study, 
two  semi-explicit  methods  are  examined:  a  second  order  semi-implicit  backwards 
difference  formula  (BDF2  SI) 

qn+1  =  \qn  -  V-1  +  \  At(2S(q)n  -  Sfa""1))  +  s\ AtL{qn+1  -  2 qn  +  q^1)  (3.26) 
o  o  o  o 

and  a  second  order  semi-implicit  leapfrog  method  (LF2  SI) 

qn+1  =  qn-L  +  2A tS{q)n  +  52AtL(vqn+1  -  qn  +  (1  -  v)qn~1).  (3.27) 


For  this  project,  a  value  of  v  =  0.6  was  chosen.  This  is  a  common  value  for 
atmospheric  models  [7].  To  illustrate  the  reasoning,  a  linear  stability  analysis  of  the 
semi-implicit  leapfrog  method  is  performed.  Again,  consider  the  equation 


dq 

~dt 


=  ikq 


with  solution 

q  =  Aeikt. 


Since  this  is  a  linear  equation,  L(q )  =  S(q),  and  the  semi-implicit  leapfrog 
method  is  reduced  to 


qn+1  =  qn~l  +  2A t[L(vqn+1)  +  L{{  1  -  v)qn~1)].  (3.28) 


Considering  the  case  where  v  =  .5,  and  again  letting  A  =  elkAt,  the  equation 
reduces  to 


15 


An+1  =  A”-1  +  At[ik\n+1  +  iky1-1) 


which  has  a  solution 


A 


±1 


'1  +  ikAt. 
1  —  ikAt 


(3.29) 


where  the  positive  value  represents  the  physical  solution  and  the  negative  value  rep¬ 
resents  the  computational  mode.  For  both  modes,  jA|  =  1,  meaning  the  solution 
is  neutrally  stable.  Again,  as  in  the  case  of  explicit  leapfrog  without  time-filtering, 
nonlinear  interactions  can  affect  the  physical  solution  and  cause  instabilities  [4]. 

Now,  consider  the  case  where  v  =  .6.  The  semi-implicit  leapfrog  method 
becomes 


qn+1  =  qn~l  +  2At[L{.6qn+1)  +  L(.4gn_1)].  (3.30) 


Again,  letting  A  =  elkAt,  the  method  is  reduced  to 


An+1  =  A”"1  +  2Af(.6zA;Ar!+1  +  AikW1-1) 


which  has  solution 


A 


± 


1  +  0.8ikAt 
1  —  1.2ikAt 


(3.31) 


Plotting  |A+|  and  |A_|,  shows  that  the  amplification  factor  has  been  damped, 
and  the  problems  associated  with  the  nonlinear  interactions  have  been  adressed. 

Still,  the  computational  mode  is  as  large  as  the  physical  solution.  Again, 
time-filtering  corrects  this  problem.  For  details,  see  Appendix  B. 
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Figure  2.  Stability  of  the  semi-implicit  leapfrog  time  integrator  with  no  time-filter. 
Figure  a)  uses  a  value  of  v  —  .5  and  Figure  b)  uses  a  value  of  v  —  .6. 
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Figure  3.  Stability  of  the  semi-implicit  leapfrog  time  integrator  with  a  time  filter 
weight  of  e  =  .05.  Figure  a)  uses  a  value  of  v  —  .5  and  Figure  b)  uses  a  value  of 
v  =  .6. 
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C.  HYBRID  EULERI AN-L AGRAN GIAN  SEMI-IMPLICIT 

TIME-INTEGRATORS 

Here,  the  HELSI  method  of  Giraldo  is  presented  [8].  For  this  method,  the 
governing  equations  are  rewritten  in  the  Lagrangian  form 


dn 

dt 

du 

dt 

d9 

dt 


- 7rV  •  U 

Cy 

—CpO'V  7T  —  gk  +  pV2u 

fiV2e. 


(3.32) 


More  simply,  this  can  be  written  as 


dq 

dt 


SL(q) 


where  Sl(q)  contains  the  source  terms  without  advection. 


The  HELSI  discretization  is 


(3.33) 


-£  =  {SL(q)-6L(q)}  +  6[L(q)].  (3.34) 

As  before,  the  terms  inside  the  braces  are  time-integrated  explicitly  while 
the  terms  inside  the  square  brackets  are  time-integrated  implicitly.  The  linearized 
operator  L  is  defined  in  the  same  way  as  before.  Applying  the  HELSI  operator  yields 


qn+1  =  E  amqn~m  +  l^t  £  (3mSL{q)n-m  +  S7AtL(  £  pmqn~m )  (3.35) 

m= 0  m= 0  m=—l 

where  the  tilde  above  q  denotes  the  quantity  is  determined  along  characteristics  in 
a  semi- Lagrangian  sense.  The  term  Hybrid  Eulerian- Lagrangian  comes  from  the  fact 
that  q  terms,  which  correspond  to  the  time  derivatives,  are  computed  in  the  La¬ 
grangian  sense,  while  the  remaining  terms  are  computed  in  the  Eulerian  sense.  The 
next  step  is  to  determine  how  to  compute  q ,  the  solution  values  along  the  character¬ 
istics. 
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The  goal  is  to  solve  the  Lagrangian  derivative 


^=0. 

dt 


(3.36) 


For  this  project,  the  Operator-Integration-Factor  Splitting  (OIFS)  method  [9] 


is  used. 


1.  Operator-Integration-Factor  Splitting  Method 

In  the  OIFS  method,  Eulcrian  substepping  is  used  to  compute  the  Lagrangian 
derivative.  The  goal  is  to  solve 


fs  =  ■  V«  (3.37) 

where  3s  =  ^,  ns  is  the  number  of  substeps,  and  with  q(s  =  0)  =  qn  to  approximate 
qn  and  q(s  =  0)  =  qn~ 1  to  approximate  qn  1 .  To  illustrate  this  procedure,  assume 
RK2  is  used  to  solve  equation  3.37.  This  yields 


~  n+ 

<7i 


1 

2 


9?+1 


9"  -  Y*”  •  V9" 

qn  -  A tun+5  ■  Vql+^ 


(3.38) 


and 


Q2 


qn~l  -  — fi”"1  •  Vqn~l 
q11-1  -  A tun~^  ■  VqT1 


(3.39) 


=  <72  ~^un-Vqn2  (3.40) 

qn2+1  =  qn2-Atun+^-Vqn2+h 
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where  qn  =  g”+1  and  qn  1  =  q^+l  and  the  velocity  held  is  defined  as  un+ 2  = 
u(x(tn+1),tn+^). 


Figure  4.  A  schematic  of  the  01FS  trajectory  computation.  The  curved  line  is  the 
actual  trajectory,  while  the  arrows  denote  the  paths  followed  to  compute  the  departure 
points  using  an  RK2  approximation.  Each  arrow  represents  a  two-stage  process. 


A  further  increase  in  the  maximum  allowable  timestep  can  be  accomplished  by 
using  higher  order  RK  methods,  or  by  substepping.  By  substepping,  each  timestep 
is  broken  up  into  smaller  timesteps,  which  are  passed  to  the  RK  method. 
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IV. 


TEST  CASES 


A.  CASE  1:  RISING  THERMAL  BUBBLE 

The  rising  thermal  bubble  is  presented  in  Giraldo  and  Restclli  [2] ,  and  is  similar 
to  the  test  proposed  by  Robert  [10]  .  The  initial  atmospheric  state  has  no  flow,  a 
constant  potential  temperature  field,  and  is  in  hydrostatic  balance  defined  by 


which  then  yields 

-  i  9 

7T  =  1 - -=Z 

CpO 

for  the  Exner  pressure.  The  following  potential  temperature  perturbation  is  then 
added  to  drive  the  motion  of  the  air 


e' 


0 

1+ cos  Of) 


for  r  >  rc 
for  r  <  rc 


where  6C  =  \  Kelvin,  nc  is  the  trigonometric  constant,  r  =  {x  —  xc )2  +  (z  —  zc )2 
with  the  following  constants:  6  =  300  Kelvin,  rc  =  250  meters,  and  (x,  z)  G  [0, 1000] 2 
meters  with  t  G  [0,600]  seconds  and  (xc,zc)  =  (500,350)  meters.  The  boundary 
conditions  for  all  four  boundaries  are  no-flux.  The  tests  are  run  on  a  10  element  x  10 
element  grid  with  10  polynomials  in  each  element.  Viscosity  is  ignored.  For  both  the 
explicit  and  semi-implicit  leapfrog  methods,  a  time-filter  weight  of  .05  is  used. 


B.  CASE  2:  LINEAR  HYDROSTATIC  MOUNTAIN  WAVE 

The  linear  hydrostatic  mountain  waves  is  originally  presented  in  Durran  and 
Klernp  [11]  and  Smith  [12].  The  atmosphere  is  isothermal,  meaning  the  temperature 
is  constant  at  T  =  250  Kelvin.  Initially,  the  atmosphere  has  a  mean  flow  of  u  =  20 
meters/second.  Imposing  hydrostatic  balance  requires 
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where  9  —  ? 

7 T 

This  gives 

__9_z 

7 r  =  e  cpt 

for  the  Exner  pressure. 

The  simulation  is  done  on  the  domain  (x,  z)  G  [0,  240000]  x  [0,  30000]  meters 
and  is  run  for  1  hour.  The  width  of  the  Agnesi  mountain  profile  is 

h(x,z)  =  - A —  (4.1) 

i  +  Oir) 

where  hc  —  1  meter,  xc  =  120,000  meters,  and  ac  =  10,000  meters.  The  Brunt- 
Vaisala  frequency  N  =  -j=.  Taking  N  =  0.0195/second  yields  >  1  which  is  in 
the  hydrostatic  range.  The  bottom  boundary  condition  is  no-flux,  while  the  remaining 
boundaries  use  a  non-reflecting  absorbing  boundary  condition,  which  works  like  a 
sponge.  A  spatial  resolution  of  20  elements  in  the  x-direction,  10  elements  in  the 
z-direction,  and  10  polynomials  per  element  is  used.  Viscosity  is  ignored.  For  both 
the  explicit  and  semi-implicit  leapfrog  methods,  a  time-filter  weight  of  .05  is  used. 

C.  CASE  3:  DENSITY  CURRENT 

The  density  current  is  proposed  in  Straka  et  al.  [13].  This  case  is  not  too  dis¬ 
similar  to  the  rising  thermal  bubble,  but  there  are  differences  in  the  size  of  the  domain 
and  the  shape  of  the  cold  cosine  bubble.  The  potential  temperature  perturbation  is 

9'  =  ^9C  [1  +  cos  (7 rcr)] 

where  9C  =  — 15  Kelvin,  7rc  is  the  trigonometric  constant,  and 

r  = 

The  simulation  is  done  on  the  domain  (x,  z )  G  [0,  25600]  x  [0,  6400]  meters  and 
is  run  for  900  seconds.  As  in  Straka  et  al.  [13],  only  half  of  the  horizonal  domain 
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is  defined.  A  no-flux  boundary  condtion  is  implemented  at  all  four  boundaries,  and 
a  dynamic  viscosity  of  p,  =  75  meters2 /second  is  used.  A  spatial  resolution  of  32 
elements  in  the  x-direction,  8  elements  in  the  z-direction,  and  8  polynomials  per 
element  is  used.  For  the  explicit  leapfrog  method,  a  time-filter  weight  of  .05  is  used. 
For  the  semi-implicit  leapfrog  method,  a  time- filter  weight  of  0.1  is  used. 


D.  CASE  4:  INERTIA-GRAVITY  WAVE 

This  nonhydrostatic  inertia-gravity  wave  is  identical  to  the  test  case  proposed 
by  Skamarock  and  Klemp  [14],  Initially,  the  atmosphere  is  uniformly  stratified,  has 
a  constant  mean  flow  of  u  —  10  meters/second  and  a  Brunt- Vaisala  frequency  of 
N  =  0.01/second.  By  definition 

JV2  =  9^(ll'») 

yields 

N2  _ 

9  =  90e~z 


where  60  =  300.  Hydrostatic  balance  is  imposed,  which  yields 

.2 


7T  =  1  + 


9 


N  , 

e  9  —  1 


cPe0N2 

Finally,  the  potential  temperature  is  perturbed  by 


9'  =  9r 


Sill 


(  TTcA 

V  hc  J 


!  +  (A? )' 


where  9C  =  0.01  Kelvin,  hc  =  10,000  meters,  ac  =  5,000  meters,  and  7rc  =  3.14159265 
is  the  trigonometric  constant.  The  simulation  is  done  on  the  domain  (x,  z )  G  [0,  300000]  x 
[0, 10000]  and  is  run  for  3000  seconds.  No-flux  boundary  conditions  are  implemented 
for  the  top  and  bottom  boundaries  and  periodic  boundary  conditions  are  used  for  the 
lateral  boundaries.  A  spatial  resolution  of  60-elements  in  the  x  direction,  2-elements 
in  the  z  direction,  and  10  polynomials  per  element  is  used.  Viscosity  is  ignored.  For 
both  the  explicit  and  semi-implicit  leapfrog  methods,  a  time-filter  weight  of  .01  is 
used. 
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V. 


RESULTS 


A.  OVERVIEW 

To  evaluate  performance  and  accuracy,  the  four  test  cases  were  run  using 
each  of  the  time-integration  methods.  Cases  1  and  4  were  run  on  a  single  processor, 
single  user  computer,  while  the  remaining  cases  were  run  on  a  dual-processor,  shared 
computer.  The  maximum  usable  timestep,  Courant  number  (C),  and  total  wallclock 
time  of  each  of  the  simulations  was  recorded.  For  each  test  case,  specific  metrics  were 
chosen  to  evaluate  accuracy.  Finally,  the  results  of  each  simulation  were  plotted  and 
visually  compared  as  another  method  of  comparison  between  the  time-integration 
methods. 

For  the  HELSI  methods,  the  OISF  method  was  performed  using  second  (RK2), 
third  (RK3),  and  fourth  (RK4)  order  Runge-Kutta  methods,  each  with  1,  2,  3,  and  4 
substeps,  for  a  total  of  12  HELSI  runs.  The  run  with  the  shortest  wallclock  time  and 
the  run  with  the  largest  Courant  number  are  used  for  the  comparison.  These  runs 
may  or  may  not  be  different,  giving  one  or  two  “best”  HELSI  simulations  for  each 
case.  When  multiple  HELSI  runs  have  the  same  Courant  number,  the  one  with  the 
shortest  total  wallclock  time  is  used.  The  results  of  all  runs  are  given  in  Appendix  C. 

B.  DETERMINING  THE  MAXIMUM  USABLE  TIME 
STEP 

For  each  of  the  time-integration  methods,  the  maximum  timestep  was  deter¬ 
mined  experimentally.  Using  linear  stability  theory  as  a  first  estimate  for  the  explicit 
methods,  an  initial  guess  was  made.  If  the  method  proved  stable  at  that  timestep, 
it  was  increased.  Otherwise  it  was  decreased.  The  processes  was  repeated  with  an 
increasingly  small  change  until  the  limit  of  stability  was  determined  to  a  high  degree 
of  accuracy.  Finally,  the  output  of  the  simulation  was  plotted  to  ensure  the  results 
were  consistent  with  what  was  expected. 
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For  the  semi-implicit  and  HELSI  methods,  a  similar  approach  was  used.  How¬ 
ever,  these  methods  have  a  tendency  to  preserve  stability,  but  not  accuracy,  for  in¬ 
creasingly  large  timesteps.  Therefore,  after  each  simulation,  the  output  was  plotted, 
and  the  numerical  results  were  compared  to  that  of  the  explicit  time-integrators. 
Using  these  comparisons,  it  was  clear  when  the  maximum  usuable  timestep  was  ex¬ 
ceeded. 


a)  b) 

Figure  5.  An  example  of  the  HELSI  time-integrator  retaining  stability  but  not  accu¬ 
racy.  In  figure  a),  case  4  is  time-integrated  with  a  timestep  of  4.75  seconds  using  the 
HELSI  method,  a  RK2  OIFS  method,  and  one  substep.  In  figure  b)  case  4  is  time 
integrated  using  a  timestep  of  5.75  seconds  using  the  same  HELSI  method.  While 
the  integration  remaines  stable,  the  result  is  obviously  unusable. 


The  benefit  of  an  increased  time  step  goes  beyond  decreasing  the  wallclock  time 
of  the  numerical  integration.  In  an  operational  NWP  model,  physical  parametriza- 
tions  are  done  after  each  time  step.  These  operations  can  account  for  as  much  as  70% 
of  the  total  computational  cost  of  the  forecast  integration  component  of  the  model 
run.  By  increasing  the  timestep,  the  total  number  of  physical  parametrization  calls 
decreases. 

For  example,  if  an  advanced  numerical  time-integration  method  can  double  the 
maximum  usable  timestep,  it  will  halve  the  number  of  computations  done  running  the 
parametrizations.  Even  if  this  method  requires  the  same  or  slightly  more  wallclock 
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time  to  do  the  time-integration  of  the  dynamics,  the  overall  effect  will  be  an  increase 
in  efficiency  due  to  the  reduction  of  computations  used  in  parametrizations. 


C.  CASE  1:  RISING  THERMAL  BUBBLE 

1.  Courant  Number  and  Computational  Cost 

a.  Explicit  Methods 

Table  I.  Case  1:  Courant  number  and  computational  cost  of  explicit  methods. 


Method 

Max  Timestep  (s) 

C 

Wallclock  Time  (s) 

BDF2 

.0018 

0.1340 

6617 

LF2 

.0052 

0.3871 

2387 

RK3 

.0097 

0.7220 

2263 

RK34 

.0120 

0.8932 

2315 

RK35 

.0148 

1.1020 

2102 

Of  the  explicit  methods,  the  SSP  RK  methods  offer  the  highest  degree  of 
efficiency  while  BDF2  underperforms  significantly.  Still,  the  Courant  number  reaches 
a  limit  around  1. 


b.  Semi-Implicit  Methods 


Table  II.  Case  1:  Courant  number  and  computational  cost  of  semi- implicit  methods. 


Method 

Max  Timestep  (s) 

C 

Wallclock  Time  (s) 

BDF2  SI 

.55 

40.94 

644.6 

LF2  SI 

.85 

63.27 

926.6 

As  expected,  the  Semi-Implicit  methods  offer  significant  improvement 
in  effeciency  over  all  of  the  explicit  methods.  LF2  SI,  for  example,  allows  for  a  Courant 
number  over  50  times  greater  than  the  best  explicit  methods.  Meanwhile,  BDF2  SI 
reduces  Wallclock  time  by  about  70%. 
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c. 


HELSI  Method 


Table  III.  Case  1:  Courant  number  and  computational  cost  of  HELSI  method. 


Depart  Method 

NS 

Max  Timestep  (s) 

C 

Wallclock  Time  (s) 

RK4 

4 

11.25 

837.4 

215.60 

For  the  HELSI  comparison,  depart  method  refers  to  the  time- integrator 
used  in  the  OIFS  method,  while  NS  is  the  number  of  substeps.  The  HELSI  method 
offers  a  great  improvement  in  efficiency  over  the  semi-implicit  methods,  with  wallclock 
time  reduced  by  a  factor  of  about  13  to  20,  and  offers  an  even  greater  improvement 
over  the  explicit  methods.  In  fact,  HELSI  completes  the  simulation  in  roughly  one 
tenth  the  time  of  RK35,  and  with  a  Courant  number  around  760  times  larger. 

2.  Accuracy 

Since  no  analytic  solution  is  available  for  this  test  case,  the  maximum  and 
minimum  values  of  the  potential  temperature  perturbations  6'  are  compared  to  ensure 
consistency  among  the  methods.  Potential  temperature  is  chosen  as  this  is  a  thermal 
problem.  Recall  that  the  maximum  temperature  perturbation  within  the  rising  bubble 
at  the  start  of  the  simulation  is  .5  degrees  Kelvin.  The  explicit  methods  are  compared 
to  one  another,  and  the  overall  most  efficient  method  is  chosen  for  comparison  with 

the  advanced  methods. 

a.  Explicit  Methods 


Table  IV.  Case  1:  A  comparison  of  potential  temperature  perturbations  for  the  ex¬ 
plicit  time-integrators. 


Method 

Qt 

umax 

Qt 

umin 

BDF2 

.513 

-.0849 

LF2 

.534 

-.0959 

RK3 

.533 

-.0959 

RK34 

.535 

-.0960 

RK35 

.535 

-.0959 

Due  to  its  extreme  computational  cost,  BDF2  is  emliminated  from  the 
comparison  immediately.  Amongst  the  other  four  methods,  the  results  are  virtually 
identical.  Therefore,  RK35  is  chosen  for  comparison  to  the  advanced  methods  because 
of  its  effeciency. 

b.  Semi-Implicit  Methods 


Table  V.  Case  1:  A  comparison  of  potential  temperature  perturbations  for  the  semi- 
implicit  time-integrators. 


Method 

nr 

umax 

nr 

umin 

BDF2 

.509 

-.0799 

LF2 

.529 

-.0519 

Both  of  the  semi-implicit  methods  give  results  in  line  with  expectations. 
BDF2  SI  is  chosen  for  comparison  as  it  has  a  shorter  wallclock  time  than  LF2  SI 

c.  HELSI  Method 


Table  VI.  Case  1:  Potential  temperature  perturbations  for  the  HELSI  time-integrator. 

Depart  Method  NS  Q'max  d'min 
RK4  4  .521  -.0911 


The  results  of  the  HELSI  simulation  are  in  line  with  those  of  the  explicit 
runs,  which  are  known  to  be  very  accurate.  The  best  methods  are  plotted  to  show 
that  the  qualitative  results  are  similar  for  each. 
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a)  b) 

Figure  6.  Casel:  A  comparison  of  potential  temperature  perturbations  using  a)  the 
RK35  explicit  and  b)  the  BDF2  semi-implicit  time-integrators  for  case  1  at  3000  sec. 


Figure  7.  Case  1:  A  view  of  potential  temperature  perturbations  using  the  HELS1 
time-integrator  with  RK4  depart  method  and  4  substeps  at  3000  sec. 

3.  Comparison  and  Conclusions 

This  is  an  excellent  example  of  the  potential  of  the  HELSI  method.  Not  only 
does  the  HELSI  method  allow  for  a  dramatically  increased  Courant  number  and  run 
significantly  faster  than  the  semi-implicit  methods,  it  also  gives  results  very  similar 
to  the  explicit  methods,  which  are  known  to  be  very  accurate. 

Figure  6  b)  shows  that  the  BDF2  SI  method  seems  to  smooth  the  underside  of 
the  thermal  bubble,  while  the  HELSI  method  in  Figure  7  seems  to  retain  more  of  the 
rigid  structure  of  the  RK35  method.  This  increase  in  accuracy  over  the  semi-implicit 
method,  coupled  with  the  dramatic  increase  in  efficiency,  makes  HELSI  the  logical 
choice  for  case  1. 
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D.  CASE  2:  LINEAR  HYDROSTATIC  MOUNTAIN  WAVE 


1.  Courant  Number  and  Computational  Cost 

a.  Explicit  Methods 


Table  VII.  Case  2:  Courant  number  and  computational  cost  of  explicit  methods. 


Method 

Max  Timestep  (s) 

C 

Wallclock  Time  (s) 

BDF2 

.07 

.0584 

2184.37 

LF2 

.18 

.150 

866.24 

RK3 

.34 

.284 

784.92 

RK34 

.43 

.359 

746.26 

RK35 

.53 

.442 

714.10 

Again,  the  SSP  RK3  methods  are  the  most  efficient  of  the  explicit 
time-integrators,  and  BDF2  is,  by  far,  the  slowest. 


b.  Semi-Implicit  Methods 


Table  VIII.  Case  2:  Courant  number  and  computational  cost  of  semi-implicit  meth¬ 
ods.  _ 


Method 

Max  Timestep  (s) 

C 

Wallclock  Time  (s) 

BDF2  SI 

6.0 

5.008 

287.59 

LF2  SI 

11.75 

9.807 

295.42 

The  semi-implicit  methods  also  behave  similarly  to  the  first  test  case. 
Here,  LF2  SI  increases  the  Conrant  number  more  than  20  fold,  while  both  methods 
run  the  simulation  in  less  than  half  the  time  of  RK35,  the  most  efficient  explicit 
method. 
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c. 


HELSI  Methods 


Table  IX.  Case  2:  Courant  number  and  computational  cost  of  HELSI  methods. 


Depart  Method 

NS 

Max  Timestep  (s) 

C 

Wallclock  Time  (s) 

RK2 

1 

13.0 

10.85 

231.87 

RK3 

1 

17.0 

14.19 

250.74 

Again,  the  HELSI  method  yields  great  results.  However,  there  appears 
to  be  a  limitation  on  the  Courant  number  which  is  not  directly  related  to  the  OIFS 
method.  Increasing  the  order  of  the  OIFS  method  and  number  of  substeps  cannot 
increase  the  Courant  number  beyond  14.19.  Therefore,  it  makes  the  most  sense  to 
use  the  simplest  OIFS  method  that  can  achieve  this  Courant  number,  as  increasing 
the  OIFS  method  further  will  only  slow  the  integration  without  benefit.  Still,  HELSI 
once  again  allows  for  a  significantly  higher  Courant  number  and  a  slighty  shorter 
wallclock  time  than  the  semi-implicit  methods. 

2.  Accuracy 

Unlike  the  other  test  cases,  an  analytical  solution  is  available  for  the  linear 
hydrostatic  mountain  wave.  This  allows  for  a  more  robust  comparison  of  the  accuracy 
of  each  of  the  time-integration  methods.  Root  mean  square  errors  of  each  of  the  four 
perturbation  variables  are  compared,  with  smaller  values  representing  better  error 
characteristics. 
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a. 


Explicit  Methods 


Table  X.  Case  2:  A  comparison  of  RMS  errors  for  the  explicit  time-integrators. 


Method 

nRMS 

u'rms 

v'rms 

O'rms 

BDF2 

4.62  x  10“7 

8.33  x  10"3 

4.64  x  10-4 

6.29  x  10“3 

LF2 

4.62  x  10“7 

8.33  x  10"3 

4.64  x  10-4 

6.30  x  10“3 

RK3 

4.62  x  10“7 

8.33  x  10"3 

4.64  x  10-4 

6.29  x  10“3 

RK34 

4.62  x  10“7 

8.33  x  10"3 

4.64  x  10-4 

6.29  x  10“3 

RK35 

4.62  x  10“7 

8.33  x  10"3 

4.64  x  10-4 

6.29  x  10“3 

Once  again  the  explicit  methods  all  yield  very  similar  results.  Since  all 
the  error  characteristics  are  virtually  identical,  RK35  is  again  chosen  for  comparison, 
as  it  is  the  most  efficient  method. 

b.  Semi-Implicit  Methods 


Table  XI.  Case  2:  A  comparison  of  RMS  errors  for  the  semi-implicit  time-integrators. 


Method 

k'rms 

u'rms 

v'rms 

/V 

URMS 

BDF2  SI 

4.61  x  10"7 

8.31  x  HT3 

4.63  x  10"4 

6.28  x  10“3 

LF2  SI 

4.61  x  10"7 

8.34  x  10"3 

4.65  x  10"4 

6.29  x  10“3 

For  this  test  case,  the  semi-implicit  methods  seem  to  perform  about  as 
well  as  the  explicit  methods. 

c.  HELSI  Method 


Table  XII.  Case  2:  A  comparison  of  RMS  errors  for  the  HELSI  time-integrators. 


Depart  Method 

NS 

Krais 

u'rms 

v'rais 

Orms 

RK2 

1 

4.54  x  lO"7 

8.27  x  10~3 

4.51  x  10-4 

6.23  x  10“3 

RK3 

1 

4.56  x  10“7 

8.50  x  10~3 

4.71  x  10“4 

6.26  x  HT3 
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The  HELSI  simulation  using  RK2  for  the  OIFS  method  actually  yields 
slightly  better  error  characteristics  than  both  the  explicit  and  semi-implicit  time- 
integrators.  The  simulation  using  RK3  for  the  OIFS  method,  however,  introduces 
slighty  larger  error  characteristics  into  the  velocity  vector.  Still,  these  errors  are  quite 
good,  and,  given  the  overall  increase  in  efficiency,  are  likely  acceptable. 


a)  b) 

Figure  8.  Case  2:  A  comparison  of  vertical  velocity  perturbations  for  a)  the  RK35 
explicit  and  b)  the  BDF2  semi-implicit  time-integrators  at  1  hr. 


Figure  9.  Case  2:  A  comparison  of  the  vertical  velocity  perturbations  for  a)  the 
HELSI  method  with  depart  method  RK2  and  1  substep  and  b)  the  HELSI  method 
with  depart  method  RK3  and  1  substep  at  1  hr. 


3.  Comparison  and  Conclusions 

Once  again  the  HELSI  time-integrator  is,  by  far,  the  most  efficient  time  inte¬ 
gration  method.  However,  this  test  case  exposes  a  limitation  of  the  HELSI  method. 
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The  Courant  number  could  not  be  increased  beyond  14.19.  The  likely  limiting  factor 
is  the  linearization  used  in  the  semi-implicit  part  of  the  HELSI  operator.  As  the 
timestep  becomes  increasingly  large,  the  linear  approximations  break  down,  and  the 
solution  cannot  converge. 

The  accuracy  of  the  HELSI  method  also  seems  to  suffer  slightly  for  this  test 
case.  While  the  RMS  errors  of  both  HELSI  simulations  are  quite  close  to  those  of  the 
explicit  runs,  Figure  9  b),  where  At  =  17  seconds,  shows  some  artifacts  on  the  right 
side  of  the  graphic  which  are  not  seen  in  the  explicit  runs.  Figure  9  a),  with  At  =  13 
seconds,  is  more  inline  with  the  RK35  method,  and  still  offers  an  improvement  in 
both  Courant  number  and  wallclock  time  over  the  semi-implicit  methods. 

Overall,  the  HELSI  method  again  seems  to  be  the  best  choice.  The  benefits 
may  not  be  as  dramatic  as  the  first  test  case,  yet  HELSI  still  offers  a  significant 
improvement. 

E.  CASE  3:  DENSITY  CURRENT 

1.  Courant  Number  and  Computational  Cost 

a.  Explicit  Methods 

Table  XIII.  Case  3:  Courant  number  and  computational  cost  of  explicit  methods. 


Method 

Max  Timestep  (s) 

C 

Wallclock  Time  (s) 

BDF2 

.02 

.1336 

1542.40 

LF2 

.0575 

.3841 

549.88 

RK3 

.115 

.7628 

607.47 

RK34 

.140 

.9352 

564.17 

RK35 

.175 

1.1690 

565.22 

For  the  density  current,  BDF2  is,  once  again,  the  least  efficient  time- 
integrator.  LF2  is  surprisingly  efficient,  with  a  wallclock  time  slightly  lower  than 
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RK35.  Still,  given  the  larger  Courant  number  of  RK35,  it  should  still  be  considered 
the  most  efficient  of  the  explicit  methods. 

b.  Semi-Implicit  Methods 


Table  XIV.  Case  3:  Courant  number  and  computational  cost  of  semi-implicit  meth¬ 
ods.  _ 


Method 

Max  Timestep  (s) 

C 

Wallclock  Time  (s) 

BDF2  SI 

.625 

4.173 

243.17 

LF2  SI 

.375 

2.504 

624.67 

LF2  SI  struggles  with  this  test  case.  The  Courant  number  is  only 
double  that  of  RK35,  while  the  wallclock  time  is  actually  higher  than  all  of  the 
explicit  methods  except  BDF2.  This  is  likely  due  to  the  inclusion  of  viscosity,  which 
is  unique  to  this  case,  as  the  leapfrog  method  is  unstable  for  diffusive  problems.  BDF2 
SI,  however,  offers  significant  improvements  over  the  explicit  methods. 

c.  HELSI  Methods 


Table  XV.  Case  3:  Courant  number  and  computational  cost  of  HELSI  methods. 


Depart  Method  NS 

Max  Timestep  (s) 

C 

Wallclock  Time  (s) 

RK3  1 

1.2 

8.012 

320.12 

Once  again  the  HELSI  methods  seem  to  have  a  limitation  independent 
of  the  OIFS  method  (see  appendix  C).  Still,  the  Courant  number  is  nearly  doubled 
over  BDF2  SI.  Wallclock  time,  however,  has  increased  by  nearly  33%. 

2.  Accuracy 

The  density  current  is  another  case  without  an  analytic  solution.  Like  the 
rising  bubble,  it  is  also  a  thermal  problem.  For  this  reason,  potential  temperature 
perturbations  are  once  again  the  choice  for  comparison. 
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a. 


Explicit  Methods 


Table  XVI.  Case  3:  A  comparison  of  potential  temperature  perturbations  for  the 
explicit  time-integrators. 


Method 

O' 

® max 

nt 

umin 

BDF2 

0.294 

-9.32 

LF2 

0.281 

-9.06 

RK3 

0.279 

-9.06 

RK34 

0.271 

-9.02 

RK35 

0.268 

-8.98 

For  the  density  current,  there  is  actually  a  bit  of  difference  in  the  results 
of  the  explicit  operators.  Since  no  analytic  solution  is  available,  it  is  impossible  to 
say  which  of  these  results  is  the  most  accurate.  Still,  it  is  expected  the  truth  should 
lie  somewhere  within  this  range. 

b.  Semi-Implicit  Methods 


Table  XVII.  Case  3:  A  comparison  of  potential  temperature  perturbations  for  the 
semi-implicit  time-integrators. 


Method 

O' 

umax 

O' 

umin 

BDF2  SI 

.289 

-9.02 

LF2  SI 

.236 

-9.15 

Of  the  two  semi- implicit  methods,  BDF2  yields  results  closest  to  those 
of  the  explicit  methods.  Given  the  difference  in  the  results  produced  by  LF2  SI,  and 
its  poor  performance  (see  table  XIV),  it  does  not  seem  to  be  a  good  option  for  this 
case. 
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c. 


HELSI  Method 


Table  XVIII.  Case  3:  Potential  temperature  perturbations  for  the  HELSI  time- 
integrator. 

Depart  Method  NS  6'max  9'min 
RK3  1  .217  -9.14 


HELSI  seems  to  be  underestimating  6'max.  Still,  as  this  is  a  cold  bubble 
problem,  9'min  is  the  most  important  value  to  consider.  Therefore,  the  results  are  still 
in  line  with  the  other  time-integrators. 


a)  b) 

Figure  10.  Case  3:  A  comparison  of  potential  temperature  perturbations  for  a)  the 
RK35  explicit  and  b)  the  BDF2  semi-implicit  time-integrators  at  900  sec. 


Figure  11.  Case  3:  Potential  temperature  perturbations  for  the  HELSI  method  with 
depart  method  RK3  and  1  substep  at  900  sec. 
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3.  Comparison  and  Conclusions 

Once  again,  the  HELSI  method  has  a  limitation  other  than  the  OIFS  method. 
For  the  density  current,  this  limitation  may  also  be  due  to  the  linearization,  or  may 
be  due  to  the  inclusion  of  viscosity.  In  the  HELSI  algorithm,  viscous  terms  are  treated 
explicitly,  thus,  they  have  a  limited  stability  region. 

The  wallclock  time  for  HELSI  is  actually  larger  than  that  of  BDF2  SI.  Still, 
HELSI  offers  a  significant  increase  in  the  Courant  number.  This  improvement  will 
probably  prove  worth  the  increase  in  wallclock  time  in  an  operational  model,  as 
parametrizations  will  only  be  called  about  half  as  often. 

Figures  10  and  11  show  no  obvious  differences  between  the  time  integration 
methods.  This  is  consistant  with  the  similar  perturbation  values  in  Tables  XVI,  XVII, 
and  XVIII.  Once  again,  the  HELSI  method  is  the  best  performer  overall,  but  the 
BDF2  SI  is  not  as  far  off  as  in  the  previous  cases. 

F.  CASE  4:  INERTIAL  GRAVITY  WAVE 

1.  Courant  Number  and  Computational  Cost 

a.  Explicit  Methods 

Table  XIX.  Case  4:  Courant  number  and  computational  cost  of  explicit  methods. 


Method 

Max  Timestep  (s) 

C 

Wallclock  Time  (s) 

BDF2 

.1 

.1573 

794.4 

LF2 

.26 

.4089 

316.6 

RK3 

.47 

.7392 

299.4 

RK34 

.58 

.9122 

299.6 

RK35 

.72 

1.1320 

298.3 

Once  again,  the  explicit  time-integrators  follow  the  same  performance 
trend.  BDF2  underperforms  signficantly,  while  RK35  is  the  most  efficient. 
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b. 


Semi- Implicit  Methods 


Table  XX.  Case  4:  Courant  number  and  computational  cost  of  semi-implicit  methods. 


Method 

Max  Timestep  (s) 

C 

Wallclock  Time  (s) 

BDF2  SI 

2.3 

3.617 

111.55 

LF2  SI 

6.3 

9.908 

168.8 

The  semi-implicit  time-integrators  once  again  offer  significant  improve¬ 
ments  over  their  explicit  counterparts.  LF2  SI  allows  for  a  Courant  number  nearly  9 
times  greater  than  that  of  RK35,  while  BDF2  SI  cuts  the  wallclock  time  in  half. 


c.  HELSI  Methods 


Table  XXI.  Case  4:  Courant  number  and  computational  cost  of  HELSI  methods. 


Depart  Method 

NS 

Max  Timestep  (s) 

C 

Wallclock  Time  (s) 

RK2 

1 

4.75 

7.470 

141.63 

RK2 

3 

9.25 

14.55 

258.47 

The  HELSI  approach,  once  again  allows  for  a  significant  increase  in 
Courant  number.  Wallclock  time,  however,  is  more  than  doubled  to  acheive  this 
increase. 

2.  Accuracy 

For  this  case,  there  is  no  obvious  choice  of  metric  to  compare  accuracy.  There¬ 
fore,  to  be  consistent,  potential  temperature  perturbations  are  used. 
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a. 


Explicit  Methods 


Table  XXII.  Case  4:  A  comparison  of  potential  temperature  perturbations  for  the 
explicit  time-integrators. 


Method 

f V 

umax 

nr 

umin 

BDF2 

.00279 

-.00152 

LF2 

.00279 

-.00152 

RK3 

.00279 

-.00152 

RK34 

.00279 

-.00152 

RK35 

.00279 

-.00152 

Since  all  of  the  explicit  methods  give  the  exact  same  potential  temper¬ 
ature  perturbations,  RK35  is  once  again  the  method  used  for  comparison,  as  it  is  the 
most  efficient. 


b.  Semi-Implicit  Methods 


Table  XXIII.  Case  4:  A  comparison  of  potential  temperature  perturbations  for  the 
semi-implicit  time-integrators. 


Method 

nr 

umax 

nr 

umin 

BDF2  SI 

.00280 

-.00152 

LF2  SI 

.00278 

-.00147 

The  semi-implicit  time-integrators  give  very  similar  results  to  the  ex¬ 
plicit  ones,  and  are  a  lot  more  efficient. 
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c.  HELSI  Methods 


Table  XXIV.  Case  4:  A  comparison  of  potential  temperature  perturbations  for  the 
HELSI  time-integrators. 


Depart  Method 

NS 

O' 

umax 

9' 

umin 

RK2 

1 

.00287 

-.00156 

RK2 

3 

.00280 

-.00152 

The  HELSI  methods  also  give  very  similar  results. 


a)  b) 

Figure  12.  Case  4:  A  comparison  of  potential  temperature  perturbations  for  a)  the 
RK35  explicit  and  b)  the  BDF2  semi-implicit  time-integrators  at  3000  sec. 
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a)  b) 

Figure  13.  Case  4:  A  comparison  of  potential  temperature  perturbations  for  the  a) 
HELSI  method  with  depart  method  RK2  and  1  substep  and  b)  the  HELSI  method 
with  depart  method  RK2  and  3  substeps  at  300  sec. 

3.  Comparison  and  Conclusions 

The  limitation  on  the  HELSI  approach  caused  by  the  semi-implicit  lineariza¬ 
tion  is  significant  for  this  case.  While  the  Courant  number  is  improved  over  the 
semi-implicit  time-integrators,  the  difference  is  not  nearly  as  great  as  in  the  other 
cases.  Meanwhile,  the  wallclock  times  for  the  HELSI  methods  are  significantly  higher 
than  those  of  the  semi-implicit  approach.  However,  in  an  operational  model,  the  in¬ 
crease  in  Courant  number  might  still  make  the  HELSI  method  more  efficient  overall 
when  parametrizations  are  included.  This  is  not  a  foregone  conclusion,  however,  as 
it  is  in  the  other  test  cases. 

Figures  12  and  13  show  that  all  the  time-integrators  give  similar  results.  Of 
the  HELSI  simulations,  the  RK2  OIFS  method  with  3  substeps  allows  for  the  largest 
Courant  number,  and,  therefore,  would  be  the  method  of  choice. 

Overall,  the  increased  Courant  number  likely  makes  HELSI  the  best  option. 
The  semi-implicit  approach,  however,  is  not  as  far  behind  as  the  previous  cases. 
Nonetheless,  both  methods  offer  dramatic  improvements  over  explicit  time-integrators. 
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VI.  CONCLUSIONS  AND 
RECOMMENDATIONS 

HELSI  time-integrators  offer  a  dramatic  improvement  in  Courant  number 
when  compared  to  semi-implicit  and  explicit  time-integrators.  Their  accuracy  ap- 
proachs  that  of  the  explicit  methods,  and  meets  or  exceeds  that  of  the  semi-implicit 
time-integrators  in  most  cases.  Considering  the  vast  improvements  in  efficiency,  the 
operational  use  of  HELSI  time-integrators  seems  to  be  a  logical  approach. 

Implementing  HELSI  time-integrators  in  an  NWP  model  would  free  significant 
computational  resources  to  be  used  in  other  areas.  They  would  allow  for  significant 
improvements  in  data  assimilation,  parametrizations,  and/or  spatial  resolution  with¬ 
out  large  increases  in  hardware  requirements.  The  overall  benefits  of  these  improve¬ 
ments  would,  in  all  likelihood,  result  in  a  better  forecast,  even  though  the  HELSI 
methods  may  be  slightly  less  accurate  for  a  given  resolution  than  currently  used 
split-explicit  methods. 

In  order  for  their  potential  to  be  fully  realized,  however,  a  few  issues  need  to  be 
addressed.  One  limitation  of  the  current  HELSI  time-integrators  is  in  the  handling  of 
viscosity.  Currently,  viscosity  is  handled  explicitly  in  the  HELSI  algorithm.  Treating 
viscosity  in  the  semi-implicit  part  of  the  algorithm  may  alleviate  some  of  the  problems 
it  causes. 

Another  possible  improvement  to  the  HELSI  approach  would  be  the  use  of  the 
RK35  scheme  in  the  OIFS  method.  Since  this  method  proved  to  be  the  most  efficient 
of  the  explicit  methods  tested,  it  is  likely  that  it  would  provide  the  best  efficiency  to 
the  OIFS  approach.  This  improvement,  however,  could  only  be  realized  if  the  other 
limitations  were  first  addressed,  as  the  HELSI  method  is  currently  limited  by  other 
factors. 

Finally,  the  large  Courant  number  of  the  HELSI  operator  exposes  a  weakness 
in  the  approach.  The  linearization  used  in  the  semi-implicit  part  of  the  HELSI  method 
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often  breaks  down  with  large  timesteps.  A  possible  remedy  would  be  to  implement  a 
fully-implicit  scheme,  which  is  not  a  trivial  task.  However,  the  potential  benefits  of  a 
fully  implicit  method  make  it  worth  continuing  to  explore  despite  the  difficulties. 

This  work  has  shown  that,  overall,  the  HELSI  approach  has  several  properties 
that  are  superior  to  those  of  current  methods,  even  with  its  shortcomings.  Addressing 
these  problems  can  only  make  HELSI  time-integrators  an  even  greater  asset  to  the 
operational  NWP  community. 
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APPENDIX  A.  COEFFICIENTS  FOR  RK35 

METHOD 


coefficient  value 

a[  0.377268915331368 
a2  0.377268915331368 
a3  0.355909775063327 
a4  0.644090224936674 
a5  0.242995220537396 
a6  0.367933791638137 
a7  0.632066208361863 
a8  0.238458932846290 
a9  0.762406163401431 
a10  0.237593836598569 
an  0.287632146308408 
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APPENDIX  B.  STABILITY  ANALYSIS  OF 
LEAPFROG  METHODS  WITH 
TIME-FILTERING 

1.  EXPLICIT  LEAPFROG  WITH  TIME-FILTER 

Starting  with  equation  3.3  and  3.4,  the  definition  of  the  time-filtered  leapfrog 
method  becomes 


qn+ 1  =  qn~l  +  2A  t{S{qn)).  (B.l) 

Consider  the  equation 

=  ikq  (B.2) 

with  solutions 

q  =  Aeikt  (B.3) 

for  the  non-filtered  solution  and 

q  =  Aeikt  (B.4) 

for  the  filtered  solution. 

If  the  equations  are  discretized,  then,  at  timestep  n,  t  =  nAt, 

qn  =  AeiknAt  =  A(eikM)n  (B.5) 


and 


qn  =  AeiknM  =  A{eikAt)n  (B.6) 

S(q)n  =  ikqn  =  ikA(elkAt)n  (B.7) 
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S{q)n  =  ikqn  =  ikA(eikAt)n. 

Consider  the  case  where  there  is  a  time-filter,  then 

A(eikAt)n+1  =  A(eikAt)n~ 1  +  2A  tikA(eikAt)n 

or 

A\n+1  =  A\n_1  +  2ikAtXn 
where  A  =  elkAi .  Dividing  by  A  yields 

AX2  =  A  +  2ikAtX 

From  the  dehnition  of  the  time-filter 

qn  =  qn  +  e(qn+1  -  2 qn  +  qn~l) 

it  is  clear  that 

AX  =  AX  +  e(AX2  -  2AX  +  A) 

or 

(A  -  e)A  =  {A  -  2Ae)X  +  AeX2. 

Multiplying  Equation  B.ll  by  (A  —  e)  yields 

A(A  —  e)A2  =  (A  —  e)A  +  2ik{X  —  e)AtX. 

Using  Equation  B.14,  this  can  be  rewritten  as 

A(A  —  e)A2  =  (A  —  2Ae)X  +  AeX2  +  2ik(X  —  e)AtX 

where  A  =  0  is  one  factor  and  can  be  removed  immediately,  as  it  is  trivial, 
dividing  by  AX  yields 


(B.8) 

(B.9) 

(B.10) 

(B.ll) 

(B.  12) 

(B.13) 

(B.14) 

(B.  15) 

(B.16) 

Therefore, 
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A2  —  eA  =  1  —  2e  +  eA  +  2ikAt\  —  2ikeAt 


(B.17) 


or 

A2  +  (— 2e  -  2ikAt)X  +  (-1  +  2e  +  2ikeAt )  (B.18) 

whose  two  solutions  are  easily  found  using  the  quadratic  equation.  The  magnitudes 
are  then  plotted  to  give  Figure  1. 

2.  SEMI-IMPLICIT  LEAPFROG  WITH  TIME-FILTER 

Starting  with  Equation  3.28 

qn+1  =  q'1-1  +  2A t[L{vqn+1)  +  L((  1  -  v)qn~1)]  (B.19) 

adding  the  time-hlter  yields 

qn+ 1  =  q’1-1  +  2 A t[L(vqn+1)  +  L((  1  -  v )qn~1)\.  (B.20) 

Using  Equations  B.5-B.8  and  the  substitution  A  =  e*fcAt,  the  time-filtered  LF2 
SI  method  becomes 

AX2  =  A  +  2At[AvikX2  +  A(  1  -  v)ik)\.  (B.21) 

Multiplying  by  (A  —  e)  yields 

AA2(A  —  e)  =  Ul(A  —  e)  +  2At[A(X  —  e)vikX2  +  Ul(A  —  e)(l  —  v)ik)\.  (B.22) 

Substituting  in  Equation  B.14.  factoring  out  AX,  and  simplifying  yields 

(1  —2ikAtv) A2  +  (—2e+AeikAtv—2ei k At) X+  (— l  +  2e+2(l  —  2e)ikAt(l—v))  (B.23) 

whose  two  solutions  can  be  found  using  the  quadratic  equation.  Plotting  their  mag¬ 
nitude  gives  Figure  3. 
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APPENDIX  C.  COMPLETE  RESULTS  FOR 

HELSI  METHODS 


Table  XXV.  Case  1:  Courant  number  and  computational  cost  of  HELSI  methods. 


Depart  Method  NS  Max  Timestep  (s)  C  Wallclock  Time  (s) 


1 

1 

74.44 

1440.8 

RK2 

2 

1.8 

134 

930.94 

3 

2.8 

208.4 

719.14 

4 

3.6 

268.0 

595.78 

1 

2.25 

167.5 

828.55 

RK3 

2 

4.00 

297.7 

499.28 

3 

6.00 

446.6 

399.15 

4 

7.75 

576.9 

352.39 

1 

3.0 

223.3 

625.37 

RK4 

2 

5.75 

428 

384.56 

3 

8.75 

651.3 

269.83 

4 

11.25 

837.4 

215.60 

Table  XXVI.  Case  2:  Courant  number  and  computational  cost  of  HELSI  methods. 


Depart  Method  NS  Max  Timestep  (s)  C  Wallclock  Time  (s) 


1 

13.0 

10.85 

231.87 

RK2 

2 

17.0 

14.19 

256.98 

3 

17.0 

14.19 

267.48 

4 

17.0 

14.19 

292.03 

1 

17.0 

14.19 

250.74 

RK3 

2 

17.0 

14.19 

277.77 

3 

17.0 

14.19 

305.28 

4 

17.0 

14.19 

328.70 

1 

17.0 

14.19 

259.56 

RK4 

2 

17.0 

14.19 

295.74 

3 

17.0 

14.19 

333.78 

4 

17.0 

14.19 

370.11 
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Table  XXVII.  Case  3:  Courant  number  and  computational  cost  of  HELSI  methods. 
Depart  Method  NS  Max  Timestep  (s)  C  Wallclock  Time  (s) 


1 

0.9 

6.01 

325.36 

RK2 

2 

1.15 

7.679 

371.15 

3 

1.2 

8.012 

398.29 

4 

1.2 

8.012 

449.39 

1 

1.2 

8.012 

320.12 

RK3 

2 

1.2 

8.012 

457.64 

3 

1.2 

8.012 

487.61 

4 

1.2 

8.012 

637.44 

1 

1.2 

8.012 

346.61 

RK4 

3 

1.2 

8.012 

479.58 

3 

1.2 

8.012 

589.83 

4 

1.2 

8.012 

701.58 

Table  XXVIII.  Case  4:  Courant  number  and  computational  cost  of  HELSI  methods. 
Depart  Method  NS  Max  Timestep  (s)  C  Wallclock  Time  (s) 


1 

4.75 

7.470 

141.63 

RK2 

2 

9 

14.15 

257.16 

3 

9.25 

14.55 

258.47 

4 

9.25 

14.55 

291.22 

1 

9 

14.15 

251.07 

RK3 

2 

9.25 

14.15 

289.94 

3 

9.25 

14.15 

326.63 

4 

9.25 

14.15 

357.31 

1 

9.25 

14.15 

271.13 

RK4 

3 

9.25 

14.15 

319.39 

3 

9.25 

14.15 

364.80 

4 

9.25 

14.15 

410.88 
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